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We  develop  a  residual-based  variational  multiscale  (RBVMS)  method  based  on  isogeo- 
metric  analysis  for  large-eddy  simulation  (LES)  of  wind-driven  shear  flow  with  Langmuir 
circulation  ( LC ).  Isogeometric  analysis  refers  to  our  use  of  NURBS  ( Non-Uniform 
Rational  B-splines)  basis  functions  which  have  been  proven  to  be  highly  accurate  in  LES 
of  turbulent  flows  (Bazilevs,  Y..  et  al.  2007,  Comput.  Methods  Appl.  Mech.  Eng..  197,  pp. 
173-201 ).  LC  consists  of  stream-wise  vortices  in  the  direction  of  the  wind  acting  as  a  sec¬ 
ondary  flow  structure  to  the  primary,  mean  component  of  the  flow  driven  by  the  wind.  LC 
results  from  surface  wave-current  interaction  and  often  occurs  within  the  upper  ocean 
mixed  layer  over  deep  water  and  in  coastal  shelf  regions  under  wind  speeds  greater  than 
3  m  s~J .  Our  LES  of  wind-driven  shallow  water  flow  with  LC  is  representative  of  a 
coastal  shelf  flow  where  LC  extends  to  the  bottom  and  interacts  with  the  sea  bed  bound¬ 
ary  layer.  The  governing  LES  equations  are  the  Craik-Leivobich  equations  (Tejada-Mar- 
tlnez,  A.  E.,  and  Grosch,  C.  E.,  2007,  J.  Fluid  Mech.,  576,  pp.  63-108;  Gargett,  A.  E., 
2004,  Science,  306,  pp.  1925-1928),  consisting  of  the  time-filtered  Navier-Stokes  equa¬ 
tions.  These  equations  possess  the  same  structure  as  the  Navier-Stokes  equations  with  an 
extra  vortex  force  term  accounting  for  wave-current  interaction  giving  rise  to  LC.  The 
RBVMS  method  with  quadratic  NURBS  is  shown  to  possess  good  convergence  character¬ 
istics  in  wind-driven  flow  with  LC.  Furthermore,  the  method  yields  LC  structures  in  good 
agreement  with  those  computed  with  the  spectral  method  in  ( Thorpe ,  S.  A.,  2004,  Annu. 
Rev.  Fluids  Mech.,  36,  pp.  584  55-79 )  and  measured  during  field  observations  in 
(D'Alessio,  S.  J.,  et  al.,  1998,  J.  Phys.  Oceanogr.,  28,  pp.  1624-1641 ;  Kantha,  L.,  and 
Clayson,  C.  A.,  2004,  Ocean  Modelling,  6,pp.  101-124).  [DOT:  10.1115/1.4005059] 


1  Introduction 

Wave  effects  play  an  important  role  in  determining  surface 
boundary  fluxes  of  momentum,  energy  and  scalars  and  ultimately 
vertical  mixing  [1].  Wave-current  interaction  is  among  several  flow 
phenomena  generating  turbulence  in  the  ocean;  others  include 
wind-and  tidal-driven  shear,  buoyancy-driven  convection  and  wave 
breaking.  Wind  speeds  greater  than  3m  s_1  often  lead  to  wave- 
current  interaction  sufficiently  strong  to  generate  Langmuir  circula¬ 
tion  (LC),  consisting  of  pairs  of  parallel  counter-rotating  vortices 
(or  cells)  oriented  approximately  in  the  downwind  direction,  as 
shown  in  the  sketch  in  Fig.  1.  The  cells  are  characteristic  of  the  tur¬ 
bulence  (i.e.,  the  Langmuir  turbulence)  advected  by  the  mean  flow. 

The  surface  convergence  of  each  cell  generates  a  downwelling 
region  characterized  by  negative  vertical  velocity  fluctuations 
while  the  bottom  divergence  generates  an  upwelling  region  charac¬ 
terized  by  positive  vertical  velocity  fluctuations,  leading  to  increased 
levels  of  mixing.  Bubbles,  particulate  matter  and  flotsam  accumulate 
along  the  surface  convergence  of  the  cells  forming  what  are  often 
referred  to  as  “windrows.”  Surface  convergences  of  the  cells  are 
characterized  by  intensification  of  positive  downwind  velocity  fluctu¬ 
ations  leading  to  an  enhanced  mean  current  as  seen  in  Fig.  1. 

Historically,  Langmuir  cells  have  been  measured  within  the 
upper  ocean  surface  mixed  layer  in  deep  water  far  above  the  bot¬ 
tom  (Fig.  1).  However,  Gargett  et  al.  (2004)  reported  detailed 
acoustic  Doppler  current  profiler  (ADCP)  measurements  of  Lang¬ 
muir  cells  engulfing  the  entire  water  column  lasting  as  long  as  18 
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hs  in  a  shallow  water  region  off  the  coast  of  New  Jersey.  Measure¬ 
ments  were  made  at  Rutgers'  LEO  15  cabled  observatory  in  15  m 
depth  water.  Gargett  et  al.  (2004)  denoted  the  observed  full-depth 
cells  as  Langmuir  supercells  (LSC)  because  of  their  important  con¬ 
tribution  towards  transport  of  sediment  and  bioactive  material  on 
shallow  shelves.  The  strong  coherence  of  LSC  makes  them  more 
effective  than  classical  bottom  boundary  layer  turbulence  at  moving 
material  out  of  the  low-speed  layer  near  bottom  and  into  the  strong 
and  strongly  directional  downwind  mean  flow  associated  with  these 
events.  Gargett  et  al.  (2004)  suggested  that  transport  in  such  super¬ 
cell  events  dominates  net  annual  transport  of  sediment  at  LE015. 

Originally  described  by  Langmuir  [2],  LC  is  now  generally 
accepted  to  be  the  result  of  wave-current  interaction  or,  more  spe¬ 
cifically,  the  interaction  between  the  wind-driven  shear  current 
and  the  Stokes  drift  current  induced  by  surface  gravity  waves  [2]. 
As  with  all  turbulence,  Langmuir  turbulence  encompasses  a  range 
of  spatial  and  temporal  scales.  Among  the  larger  spatial  scales  are 
those  of  the  cells  which  extend  in  the  downwind  direction  for  tens 
of  meters  to  kilometers  and  are  separated  by  distances  ranging 
from  meters  to  as  much  as  a  kilometer  [3], 

A  mechanism  for  the  generation  of  LC  was  first  proposed  by 
Craik  and  Leibovich  [4].  It  consists  of  a  vortex  force  (the  Craik- 
Leibovich  force  or  C-L  force)  in  the  momentum  equations  repre¬ 
senting  the  interaction  between  the  Stokes  drift  velocity,  induced 
by  surface  gravity  waves,  and  the  vertical  shear  of  the  current; 
specifically,  the  C-L  vortex  force  is  the  vector  cross  product 
between  the  Stokes  drift  velocity  and  the  vorticity  of  the  flow. 
Main  parameter  ingredients  in  this  force  are  the  dominant  wave¬ 
length  and  amplitude  of  the  surface  waves  used  to  define  the 
Stokes  drift  velocity  profile.  The  square  root  of  the  ratio  of  wind 
friction  velocity  to  surface  Stokes  drift  velocity  defines  the 
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Fig.  1  Sketch  of  Langmuir  circulation 

turbulent  Langmuir  number,  La, ,  a  measure  of  the  strength  of 
wind-driven  shear  forcing  relative  wave  forcing  [5],  The  strength 
of  the  Langmuir  cells  is  thus  inversely  proportional  to  La,. 
Although  this  definition  is  counterintuitive,  it  is  widely  used  within 
the  LES  community  [6].  A  more  intuitive  definition  is  described  in 
Ref.  [7]  as  the  ratio  of  Stokes  production  of  turbulence  kinetic 
energy  (TKE)  to  shear  production  of  TKE.  In  this  definition,  the 
strength  of  the  Langmuir  turbulence  increases  with  La,. 

The  C-L  force  arises  from  low-pass  time  filtering  or  wave-phase 
averaging  of  the  Navier-Stokes  equations  in  order  to  filter  out  the 
high  frequency  surface  waves.  Hereafter,  the  time  filtered  Navier- 
Stokes  equations  with  the  C-L  force  will  be  referred  to  as  the  C-L 
equations.  Inclusion  of  the  C-L  force  in  the  momentum  equation 
greatly  reduces  the  computational  complexity  as  it  eliminates  the 
need  to  resolve  the  surface  waves  giving  rise  to  LC.  Instead,  the  top 
of  the  flow  domain  is  simply  taken  to  be  bounded  by  a  flat  (nonde¬ 
forming)  surface  denoting  the  mean  water  height.  Note  that  the  C-L 
framework  does  not  account  for  the  impact  of  wave-breaking  on  the 
turbulence  resolved.  Sullivan  and  McWilliams  [8]  have  incorporated 
a  stochastic  model  of  wave  breaking  to  the  C-L  equations  in  their 
LES  of  Langmuir  circulations  within  the  upper  ocean  mixed  layer. 

The  C-L  equations  have  enabled  a  number  of  successful  LES 
(large-eddy  simulations)  and  RANS  (simulations  based  on  the 
Reynolds  averaged  Navier-Stokes  equation)  describing  the  verti¬ 
cal  and  horizontal  structure  of  upper  ocean  Langmuir  turbulence 
in  statistical  equilibrium,  e.g.,  [5,9,10-13].  See  the  review  in  Ref. 
[10]  for  further  references.  In  direct  numerical  simulation  (DNS) 
all  of  the  scales  of  a  turbulent  flow  are  explicitly  computed  while 
in  RANS,  only  the  mean  component  of  the  flow  is  explicitly  com¬ 
puted  and  the  turbulent  scales  or  eddies  are  parameterized.  Large- 
eddy  simulations  is  a  compromise  between  the  computationally 
expensive  DNS  and  the  less  demanding  RANS.  In  LES  the  more 
energetic,  larger  eddies  are  explicitly  computed  while  the  smaller 
eddies  are  modeled  [14],  Note  that  a  DNS  of  Langmuir  turbulence 
would  not  be  possible  with  the  C-L  equations,  as  these  equations 
are  the  result  of  time  filtering.  A  DNS  would  require  a  prohibit- 
edly  expensive  simulation  resolving  surface  waves  and  thus  direct 
resolution  of  the  interaction  between  surface  waves  and  the  shear 
current  giving  rise  to  Langmuir  circulation. 

Of  particular  interest  here  is  the  work  of  Tejada-Martinez  and 
Grosch  [9]  who  performed  LES  of  a  finite-depth,  wind-driven 
shear  current  with  full-depth  Langmuir  circulation  guided  by  the 
measurements  in  Refs.  [15,16].  Their  LES  used  the  C-L  equations. 
In  particular,  predictions  from  the  LES  compared  favorably  with 
in-water  measurements  of  the  Reynolds  stress  components.  Fur- 
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thermore,  the  structure  of  the  Langmuir  cells,  manifested  as  a  sec¬ 
ondary  turbulent  structure  in  the  shear  flow,  was  consistent  with 
that  measured  during  the  in-water  observations. 

Ideally,  computational  methods  for  turbulent  flow  must  be  able 
to  resolve,  or  accurately  model  all  the  relevant  flow  scales  and  their 
interactions  in  the  presence  of  complex  geometrical  configurations. 
Currently,  computational  approaches  for  turbulent  flows  include  ef¬ 
ficient  techniques  that  are  based  on  high-order  functions  (e.g., 
global  polynomial  and  Fourier  bases)  able  to  accurately  represent 
the  small  scale  physics  in  turbulent  processes  at  relatively  high 
Reynolds  numbers,  albeit  on  very  simple  geometrical  configura¬ 
tions.  For  example,  LES  performed  with  the  C-L  equations  have 
focused  on  simple  flow  configurations  while  accurately  resolving 
complex  physical  phenomena  such  as  upper  ocean  mixed  layer 
(UOML)  turbulence  and  internal  waves  in  the  pycnocline  below 
the  UOML  [17,18].  On  the  other  end  of  the  spectrum,  much  pro¬ 
gress  has  been  made  over  the  last  several  decades  on  the  computa¬ 
tion  of  flows  over  complex  geometrical  configurations.  These 
methods,  based  on  low-order  functional  representations,  are  able  to 
represent  relatively  well  the  gross  features  of  a  given  turbulent 
flow,  yet  they  do  not  possess  the  high-order  accuracy  of  the  afore¬ 
mentioned  spectral  techniques  to  predict  the  more  detailed  features 
of  the  flow.  Thus,  it  appears  that  there  is  a  gap  between  techniques 
capable  of  accurately  resolving  all  scales  in  turbulent  flows  on  sim¬ 
ple  geometries  and  techniques  capable  of  resolving  only  large  scale 
(gross)  features  on  complex  geometries.  In  order  to  bridge  this  gap, 
a  methodology  is  necessary  that  simultaneously  possesses  superior 
approximation  and  uniform  convergence  behavior  over  a  wide 
range  of  spatial  and  temporal  scales,  necessary  for  capturing  flow 
physics,  and  the  geometrical  flexibility,  necessary  for  geophysical 
and  engineering  applications.  We  feel  that  isogeometric  analysis 
based  on  Non-Uniform  Rational  B-splines  (NURBS),  recently  pro¬ 
posed  as  a  new  computational  technique  by  Hughes  et  al.  [19],  is 
an  excellent  candidate  for  the  task.  NURBS  are  spline  basis  func¬ 
tions  used  for  representation  of  complex  geometry,  are  locally 
supported,  and  possess  spectrallike  approximation  properties  com¬ 
pared  to  standard  complex-geometry  approaches  (i.e.,  low-order  fi¬ 
nite  elements  and  finite  volumes). 

In  this  paper,  we  extend  variational  multiscale  turbulence  model¬ 
ing  procedures  to  the  C-L  equations  and  implement  them  using 
isogeometric  analysis  based  on  NURBS.  In  the  case  of  the  Navier- 
Stokes  equations,  this  framework  has  been  shown  to  yield  spectral¬ 
like  resolution  while  rapidly  converging  under  refinement  to  direct 
numerical  simulation  (DNS),  and  providing  good  quality  large-eddy 
simulation  (LES)  results  on  intermediate  meshes  (e.g.,  see  Bazilevs 
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et  al.  [20]),  for  both  homogeneous  as  well  as  boundary  layer  and  tran¬ 
sitional  flows.  Extension  of  these  methodologies  to  the  C-L  equations 
will  allow  future  research  investigating  effects  of  complex  lateral 
boundaries  and  bathymetry  on  Langmuir  circulation  and  other  salt 
and  momentum  mixing  processes  in  coastal  oceans  and  estuaries. 

Solution  of  the  Navier-Stokes  equation  augmented  with  the  C-L 
vortex  force  is  nontrivial  as  the  latter  term  is  an  advective  term 
giving  rise  to  instabilities  requiring  stabilization  of  the  type 
presented  in  this  paper.  Note  that  the  C-L  vortex  force  has  been 
previously  identified  in  Refs.  [10]  and  [21]  as  giving  rise  to  insta¬ 
bilities  by  triggering  scales  of  size  smaller  than  the  grid  (i.e.,  sub- 
grid-scales).  The  goal  of  this  paper  is  the  nontrivial  extension  of 
advection  stabilization  within  a  variational  formulation  in  order  to 
consistently  account  for  the  advective  nature  of  the  C-L  vortex 
force.  Advection  stabilization  behaves  as  a  subgrid-scale  model 
derived  from  the  Navier-Stokes  equations  (in  our  case  the  Navier- 
Stokes  equations  augmented  with  the  C-L  vortex  force)  via  a  split¬ 
ting  of  the  space  of  solutions.  The  splitting  results  in  finite¬ 
dimensional  equations  governing  the  large-eddies  and  infinite 
dimensional  equations  governing  unresolved  eddies.  Analytical 
expressions  for  the  solution  of  the  latter  are  obtained  via  mathe¬ 
matical  approximations  [1],  The  splitting  of  the  space  of  solutions 
does  not  involve  application  of  spatial  filtering  nor  does  it  gener¬ 
ate  a  residual  subgrid-scale  stress,  as  is  the  case  in  classical  LES. 
In  this  approach  the  subgrid-scale  model  is  not  for  the  subgrid- 
scale  stress,  but  rather  for  a  subgrid-scale  velocity.  Thus,  a  con¬ 
ventional  subgrid-scale  stress  is  not  present  in  the  formulation.  As 
noted  earlier,  Bazilevs  et  al.  [20]  have  shown  that  this  approach 
yields  LES  solutions  on  intermediate  meshes  while  rapidly  con¬ 
verging  to  direct  numerical  simulation  (DNS)  solutions  on  finer 
meshes.  Their  simulations  involved  turbulent  channel  flows  and 
forced  isotropic  turbulence. 
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and  summation  on  the  repeated  index  i  is  employed.  In  the  follow¬ 
ing  section  we  present  the  numerical  formulation  of  the  above  par¬ 
tial  differential  equations. 


2.2  The  Space-Discrete  Formulation  of  the  Craik- 
Leibovich  Equations.  In  this  section  we  present  the  residual- 
based  variational  multiscale  (RBVMS)  formulation  of  the  C-L  par¬ 
tial  differential  equations.  The  formulation  is  a  straight-forward 
extension  of  the  RBVMS  formulation  for  the  incompressible 
Navier-Stokes  equations  [20,22]  that  also  accounts  for  the  presence 
of  the  C-L  forcing  terms.  For  better  approximation  of  thin  bound¬ 
ary  layers  near  no-slip  walls,  weak  enforcement  of  the  Dirichlet 
boundary  conditions,  proposed  in  Ref.  [23],  is  also  employed. 

Let  V'  denote  the  discrete  solution  space  for  the  velocity- 
pressure  pair  { uh,  ph  j  and  let  Wh  denote  the  discrete  weighting 
space  for  the  linear  momentum  and  continuity  weighting  functions 
\du\  dph). 

The  space-discrete  Navier-Stokes  problem  is  stated  as:  Find 
{ n\  ph  ]  £  V1'  such  that  5p7'}  £  VV*, 


2  Numerical  Approach 

In  the  following  subsections  we  provide  precise  details  of  the  nu¬ 
merical  method  deployed  for  solving  the  Craik-Leibovich  equations. 

2.1  The  Craik-Leibovich  Equations  at  the  Continuous 
Level.  Let  fi  £  K3  be  the  problem  domain  and  let  d £2  denote  its 
boundary.  A  conservative  form  of  the  dimensionless  Navier- 
Stokes  equations,  with  C-L  forcing,  in  the  Eulerian  frame  are 
taken  as  a  starting  point  of  our  developments,  and  are  given  as 

^  +  V  ■  (u  ®  «)  +  Vp  —  V  ■  (2 Re^'Vsu) 

—  La;2cj>  x  V  x  u  =  0  in  Q  (1) 

V  ■  u  =  0  in  Q  (2) 


B({Su\  Sph } ,  { uh , ph } )  +  Bvms ({ duh ,  Sph } ,  { uh ,ph}) 

+  Bwbc({Suh,dph},{uh,ph})  =  (Su'\f)n  +  (<5w,  h)rA  (7) 

In  the  above,  (•,  -)A  denotes  an  L2-inner  product  over  A  and  F /,  is  a 
part  of  the  domain  boundary  where  traction  h  is  applied.  In  this 
work,  the  last  term  on  the  right-hand-side  of  Eq.  (7)  is  active  and 
models  the  effect  of  wind  stress  on  the  water  surface.  The  rest  of 
the  terms  of  the  above  formulation  are  defined  in  what  follows. 

B{b v-q}-{u,p})  =  (Vw,u®u)n+(q,\7  -u)n 

-(V  ■»./>)«+(".  A,-g)Q 

+  (Vsw,2Re-lVsu)n  (8) 


Equations  (1)  and  (2)  represent  conservation  of  linear  momentum 
and  mass,  respectively,  assuming  density  is  constant.  In  the  above 
u  and  p  are  the  fluid  velocity  and  pressure  (divided  by  density), 
V!  =  (1/2)(V  +  (V)r)  is  the  symmetric  spatial  gradient  of  the  ve¬ 
locity,  0  is  the  Stokes  drift  velocity  and  f  is  the  body  force  per  unit 
mass.  The  physical  nature  of  the  flow  is  characterized  by  the 
dimensionless  Reynolds  and  turbulent  Langmuir  numbers,  Re  and 
La„  respectively.  The  turbulent  Langmuir  number  is  inversely  pro¬ 
portional  to  the  strength  of  Langmuir  circulation  (i.e.,  wave  forc¬ 
ing)  relative  to  the  shear  flow  (i.e.,  wind  forcing)  in  the  model. 

The  last  term  on  the  left-hand-side  of  Eq.  (1)  represents  C-L 
forcing.  Because  the  term  depends  on  the  first-order  derivatives  of 
the  velocity  field,  it  has  the  mathematical  structure  of  advection. 
With  this  in  mind,  we  re-write  the  C-L  momentum  equations  as 

+  V  ■  (u  0  u)  +  Vp  —  V  •  (2Re  x  Vsu)  +  A =  f  in  £2 

C/t  UXj 

(3) 

where  A,  ’s  are 


is  the  Galerkin  part  of  the  weak  form 

<?},{«, p})  =  -(Vw,  «' ®  U  +  u  ®  «' +  «' 0  «% 

(9) 

are  the  RBVMS  terms,  and  the  pair  { v! ,  p'\  denotes  the  velocity 
and  pressure  subgrid  scales  (i.e.,  the  scales  that  are  too  small  to  be 
reasonably  approximated  on  a  given  mesh). 

Analogously  to  Ref.  [20],  the  subgrid  scales  are  modeled  as 

,  fdu  -  du  \ 

u  =  —  tm  [ +  Vp  —  Re  Am  +  A,- - f 

\  at  ax,  J 

p'  =  — r cV  ■  u  (10) 

where  tm  and  tc  ^  the  subgrid-scale  parameters  defined  in  what 
follows.  The  subgrid-scale  parameters  are  also  known  as  the  stabi¬ 
lization  parameters  due  to  the  similarities  between  RBVMS  and 
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stabilized  methods  [24-26].  In  fact,  stabilized  methods  are  the 
progenitors  of  the  RBVMS  methodology.  Note  that  the  fine  scales 
are  proportional  to  the  residual  of  the  C-L  partial  differential 
equations,  which  renders  the  numerical  methodology  consistent. 
Also  note  that  the  momentum  equation  residual  is  in  the  advective 
form,  which  may  be  derived  from  the  conservative  form  given  by 
Eq.  (1)  using  the  incompressibility  constraint  from  Eq.  (2). 

To  define  the  subgrid-scale  parameters  for  the  C-L  equations, 
the  momentum  equations  are  written  in  the  form  of  a  quasi-linear 
advective-diffusive  system  as 

^  +  A,  ^  —  RelAu  »f  —  Vp  in  £2,  (11) 

Ot  OX, 

where  A,  =  (u, -I  +  A,)  are  the  advective  flux  jacobians  given  by 


robustness  of  the  numerical  method  for  small  time  step  sizes  [34]. 
However,  we  did  not  pursue  this  approach  here. 

Finally,  Bwhc,  defined  as 

Bwhc{{w,q},  {u,p})  =  (w,  -2 Re-lX7su  ■  n)Fj 

+  (— 2Re~' W  ■  n,  («  —  g))Fg 

+  (w,rB(w-g))Fj  (19) 

contains  terms  that  weakly  impose  the  boundary  condition  u  =  g. 
In  Eq.  (19),  Tg  is  the  Dirichlet  part  of  the  problem  domain  bound¬ 
ary  and  g  is  the  prescribed  flow  velocity  vector.  In  the  formulation 
we  assume  that  the  normal  component  of  the  flow  velocity  vector 
is  imposed  strongly.  To  ensure  numerical  stability  and  optimal 
convergence,  the  penalty  parameter  tb  in  Eq.  (19)  is  chosen  as 


/  u  -02  -03  \ 

Ai  =  0  u  +  0!  0  I 

\0  0  u+(j)x) 

/  V  +  02  0  0  \ 

A2  =  -01  v  -03 

V  o  o  v + 02  y 


xB  =  ChRe  1  yJriiGijrij  (20) 

(12) 

where  n,’s  are  the  Cartesian  component  of  the  unit  outward  nor¬ 
mal  vector  to  rg  and  Ch  is  an  elementwise  constant  emanating 
from  a  boundary  inverse  estimate  [23],  Further  discussion  and 
computational  results  employing  weakly-enforced  Dirichlet  con¬ 
ditions  may  be  found  in  Refs.  [23,32,35,36]. 


/  w  +  03  0  0  \ 

A3  =  0  w  T  03  0  .  (14) 

\  -01  “02  w) 

The  C-L  forcing  contributions  render  these  jacobians  nondiagonal 
and  nonsymmetric,  which  requires  an  appropriate  definition  of  the 
subgrid-scale  parameters.  Based  on  the  developments  in  Ref. 
[27-29]  for  advective-diffusive  systems,  they  may  be  computed 
as  follows: 

™  =  (^2 1  +  GijAjAj  +  CtCyG.jRe  2lj  '  (15) 


where 


ij  dxj  dxj 


(16) 


is  the  metric -tensor  of  the  mapping  from  the  parametric  to  the 
physical  domain  of  the  NURBS  element,  and  Cj  is  a  constant  aris¬ 
ing  in  the  element-level  inverse  estimate  (see,  e.g.,  Ref.  [30]). 

The  definition  of  rM  in  Eq.  (15)  necessitates  the  computation  of 
the  matrix  square  root  inverse.  We  do  so  using  the  Denman- 
Beavers  algorithm  [31],  which  computes  the  matrix  square  root 
inverse  in  an  iterative  fashion.  The  algorithm  is  started  by  setting 
Xo  =  Tm  2  and  Y0  =  I,  and  the  iteration  consists  of  the  following 
updates: 


X*+i  =  x  (Xa  +  Y^T1) 

I  (17) 

Y*+t  =  2  (Y*  +  XjrT1) 

where  k  is  the  iteration  index.  In  a  small  number  of  iterations  Y 
converges  to  xM  defined  by  Eq.  (15).  Given  rM,  zc  is  computed  as 

Tc  =  (GijTMij)-'  (18) 

which  is  a  generalization  of  the  relationship  given  in  Refs.  [20,32] 
An  alternative  definition  of  the  subgrid-scale  parameters  rM 
and  tc  based  on  the  element-level  matrices  and  vectors  may  be 
employed  in  our  discrete  formulation.  The  subgrid-scale  parame¬ 
ter  definitions  given  in  Ref.  [33]  may  be  naturally  extended  to  the 
C-L  equations,  which  is  likely  to  improve  the  accuracy  and 


3  Computational  Setup 

The  computational  domain,  depicted  in  Fig.  2,  is  a  rectangular 
box  with  dimensions  2n5x2Sx  (8 /3)ti<5,  in  the  stream-wise  or 
downwind  (jc),  wall-normal  or  vertical  (y),  and  span-wise  or  cross- 
wind  (z)  directions,  respectively.  The  half-depth  of  the  domain  (in 
the  y-direction)  is  <5.  The  flow  is  driven  by  a  wind  stress  in  the  jc 
direction  applied  at  the  top  surface  (y  =  2<5),  where  no-penetration 
boundary  condition  is  also  assumed  to  hold.  No-slip  conditions 
(enforced  weakly)  are  applied  at  the  bottom  wall  boundary  (y  =  0). 
In  the  stream-wise  and  span-wise  directions  periodic  boundary 
conditions  are  employed  in  order  to  represent  an  unbounded  do¬ 
main  in  these  directions,  approximating  a  continental  shelf  flow 
miles  away  from  coastal  boundaries. 

Characteristic  flow  velocity,  Stokes  drift  velocity  and  length  are 
taken  as  wind  stress  friction  velocity  uz,  Stokes  drift  velocity  at 
the  surface  us,  and  half-depth  5,  respectively.  Characteristic  time 
scale  is  taken  as  5/uz.  Using  these  scales  to  nondimensionalize  the 
C-L  equations  gives  rise  to  the  Reynolds  number  defined  as 
Re  =  uzS/v  (where  v  is  kinematic  viscosity)  and  the  turbulent 
Langmuir  number  defined  as 

La,  =  uz/us 

The  nondimensional  Stokes  drift  0  in  Eq.  (1)  is  given  by 


(  cosh(2)cy)  ) 
2  sinh(2K<5) 

0 

V  o  / 


ye  [0,25] 


(21) 
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where  ic  =  2n/y  is  the  dominant  wave  number  and  y  is  the  domi¬ 
nant  wavelength  of  surface  gravity  waves  generating  Langmuir 
circulation  (see  Refs.  [37,38]  for  details).  In  their  LES  of  Lang¬ 
muir  circulation  within  the  upper  ocean  mixed  layer  in  deep 
water,  Harcourt  and  DAsaro  [39]  have  used  a  Stokes  drift  velocity 
in  the  C-L  force  representative  of  wind  seas  characterized  by 
broadband  equilibrium  displacement  that  represent  the  cumulative 
effect  over  some  fetch  or  duration  of  surface  stress  due  to  local 
wind  conditions.  Their  broadband  wave  spectra  have  Stokes  drift 
composed  of  a  Stokes  drift  spectrum  with  elements  that  decay 
exponentially  below  the  surface. 

The  flow  is  driven  purely  by  a  wind  stress;  thus  the  body  force, 
f,  in  Eq.  (1)  is  0.  The  dimensionless  wind  stress,  h.  is  given  by 


We  note  that  setting  <j>m  0  and  keeping  the  same  boundary  condi¬ 
tions  gives  a  similar  flow  configuration  similar  to  traditional  Cou- 
ette  flow,  but  with  a  surface  wind-stress  and  a  stationary  no-slip 
bottom  instead  of  the  usual  no-slip  bottom  and  top  plates  moving 
in  opposite  direction  to  each  other. 

For  the  computations  presented  here,  we  set  Re  =  395,  La, 
=  0.7,  and  X  =  12<5.  These  are  identical  to  the  parameters  reported 
in  spectral  LES  calculations  of  Ref.  [9],  Furthermore,  parameters 
La,  =  0.1  and  X  =125  are  characteristic  of  the  wind  and  wave 
forcing  conditions  during  the  field  measurements  of  shallow 
water,  full-depth  Langmuir  circulation  of  Gargett  and  Wells  [16]. 
Their  measurements  were  made  in  a  15  ms-deep  water  column  on 
the  coastal  shelf  off  southern  New  Jersey. 

Quadratic  NURBS  that  are  C’-continuous  across  mesh  knots 
are  employed  in  the  computations.  We  perform  our  simulations 
using  a  sequence  of  /(-refined  meshes  to  ensure  convergence  of 
the  computational  results.  The  coarsest  mesh  is  comprised  of 
24  x  26  x  24  NURBS  elements,  while  the  finest  mesh  has 
64  x  66  x  64  NURBS  elements.  In  general,  for  NURBS  of  order  p 
and  maximal  continuity  p  —  1,  the  number  of  basis  functions  in 
each  tensor-product  direction  equals  to  n+p,  where  n  is  the  num¬ 
ber  of  elements  in  this  direction.  (For  periodic  boundary  condi¬ 
tions,  the  number  of  basis  functions  is  n,  which  is  independent  of 
the  polynomial  order.)  This  is  in  contrast  to  the  C°-continuous  fi¬ 
nite  elements  of  order  p,  where  the  number  of  basis  functions  is 
pn  +  1  (or  pn  in  the  periodic  case).  This  amounts  to  significant 
savings  in  the  number  of  degrees  of  freedom  for  NURBS  elements 
with  respect  to  finite  elements  of  the  same  order,  especially  in  3D. 

The  mesh  is  uniform  in  the  periodic  directions.  The  elements  in 
the  wall-normal  (vertical)  direction  are  stretched  toward  the  top 
and  bottom  boundaries.  The  mesh  knots  are  placed  according  to 

y,  =  /u'tanh^^-^— -  -  lj  tanh_1(h)j  i  6  [l,JVy  +  1]  (23) 
where  b  =  0.973  is  used 

The  flow  is  advanced  in  time  using  the  generalized-a  method 
[40,41]  and  the  resulting  nonlinear  equations  are  solved  using  an 
inexact  Newton  Krylov  approach  (see  Ref.  [20]  for  details). 
Details  of  the  mesh  and  time  step  sizes  may  be  found  in  Table  1. 

4  Numerical  Results 

In  the  following  sub-sections  we  present  results  from  LES  of 
wind-driven  flow  with  and  without  Langmuir  circulation  (LC).  In 
the  description  of  these  results,  components  of  the  velocity  vector  u 
are  referred  to  as  either  (tq,  u2,  u3)  or  ( u .  v,  tv)  where  tq  ( u ),  u2  (v) 
and  u3  (w)  are  velocity  components  in  the  stream-wise  (downwind), 
wall-normal  (vertical)  and  span-wise  (crosswind)  directions, 
respectively. 


Table  1  Summary  of  mesh  and  time  step  sizes  used  in  the  sim¬ 
ulations.  In  the  table,  Nx,  Ny,  and  Nz  are  the  number  of  basis 
functions  used  in  the  simulation  in  each  tensor  product  direc¬ 
tion  and  Ntot  is  their  total  number.  jq+  Is  the  size  of  the  first  ele¬ 
ment  in  the  wall-normal  direction  In  non-dimensional  wall  units 
(jq+  =  utAy/v).  Time  step  size  Af  has  been  made  dimensionless 
with  characteristic  time  scale  given  as  Slu z. 


Nx 

Ny 

IV: 

N,0, 

yt 

At 

24 

26 

24 

14976 

4.62 

0.025 

32 

34 

32 

34816 

3.31 

0.0188 

48 

50 

48 

115200 

2.11 

0.0125 

64 

66 

64 

270336 

1.55 

0.00935 

4.1  Flow  Structures.  Figure  3  shows  an  instantaneous  three- 
dimensional  snapshot  of  iso-contours  of  vertical  velocity  fluctua¬ 
tions  in  the  wind-driven  flow  with  LC.  Vertical  velocity  fluctua¬ 
tions  are  characterized  by  negative  and  positive  downwind 
elongated  regions,  corresponding  to  the  downwelling  and  upwell- 
ing  limbs  of  the  Langmuir  cells  sketched  in  Fig.  1 .  Figure  4  shows 
an  instantaneous  snapshot  of  downwind  velocity  fluctuations  in 
wind-driven  flows  with  and  without  LC.  In  both  flows,  downwind 
velocity  fluctuation  is  characterized  by  downwind  elongated 
streaks  alternating  in  sign  in  the  crosswind  direction.  The  simula¬ 
tion  with  LC  was  initialized  by  “turning  on”  the  C-L  vortex  force 
in  the  simulation  without  LC  after  the  latter  had  achieved  statisti¬ 
cal  equilibrium.  Animations  (not  shown)  reveal  that  the  vortex 
force  causes  the  positive  streaks  in  the  flow  without  LC  to  merge 
together  leading  to  a  single  pair  of  streaks  (positive  and  negative). 
The  crosswind  extent  of  the  resulting  positive  streak  is  greater 
than  the  negative  streak. 

In  order  to  reveal  the  crosswind-vertical  structure  of  the  previ¬ 
ously  described  downwind  elongated  streaks,  we  perform  the  fol¬ 
lowing  triple  decomposition  of  the  computed  velocity: 

ui  =  («/)  +  {u[)lx+u"  (24) 


where  {-)a  denotes  averaging  in  time  and  over  the  downwind  ( x ) 
direction  and  the  instantaneous  velocity  fluctuation  is  obtained  via 
the  classical  Reynolds  decomposition: 

ut  =  (; Uj }  +  u\  (25) 

where  brackets  without  subscripts  denote  averaging  over  down¬ 
wind  (x)  and  crosswind  (z)  directions  and  over  time.  Note  that  the 
prime  notation  used  here  to  define  the  resolved  turbulence  veloc¬ 
ity  fluctuation  is  different  from  the  prime  notation  in  Eqs.  (9)  and 
(10)  used  to  define  unresolved  (subgrid)  scales.  The  middle  term 
on  the  right  hand  side  of  Eq.  (24)  is  defined  as  a  partially  averaged 
fluctuation: 

v'i(y,z)  =  (26) 

This  partially  averaged  velocity  fluctuation  emphasizes  coherent, 
secondary  flow  structures  in  the  downwind  direction  such  as  the 
downwind  elongated  streaks  observed  in  Fig.  4.  Figures  5  and  6 
show  the  crosswind-vertical  structure  of  the  partially  averaged  ve¬ 
locity  fluctuation  in  the  flows  with  and  without  LC,  respectively. 
Overall,  both  cases  exhibit  positive  and  negative  crosswind  cell 
structures  in  each  of  the  partially  averaged  fluctuating  velocity 
components;  the  flow  with  LC  has  a  spanwise  one-cell  structure 
while  the  flow  without  LC  has  a  less  coherent  spanwise  two-cell 
structure. 

The  one-cell  structure  in  the  flow  with  C-L  forcing  (Fig.  5)  is 
nearly  identical  to  that  obtained  with  the  spectral  LES  of  Tejada- 
Martfnez  and  Grosch  [9]  with  the  same  wind  and  wave  forcing 
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Fig.  3  Instantaneous  snapshot  of  iso-contours  of  wall-normal  (vertical)  velocity  fluctuations  in 
flow  with  LC  on  the  64  x  66  x  64  quadratic  NURBS  mesh  described  earlier 


parameters  described  earlier.  Recall  that  these  parameters  in  the 
C-L  vortex  force  have  been  chosen  as  La,  =  0.7  and  y  =  12<5  follow¬ 
ing  the  field  measurements  of  Ref.  [16].  Additionally,  the  one-cell 
structure  in  the  flow  with  C-L  forcing  possesses  all  of  the  basic 
characteristics  of  full-depth  Langmuir  circulation  expected  based 


on  the  field  measurements  in  Ref.  [16].  Experimental  data  in  Ref. 
[16]  shows  that  the  spanwise  (crosswind)  length  of  one  Langmuir 
cell  is  in  the  range  of  6 S  and  12<5.  Accordingly,  our  computation 
has  predicted  the  generation  of  only  one  Langmuir  cell  as  expected, 
given  the  crosswind  extent  chosen  for  the  domain  (see  Fig.  2). 


Flow  without  LC 


0  1  2  3  4  5  6 


Flow  with  LC 


Fig.  4  Color  maps  of  instantaneous  downwind  velocity  fluctuation  ur  on  the  downwind- 
crosswind  plane  at  mid-depth  in  flows  with  and  without  C-L  vortex  forcing  (i.e.,  with  and  without 
LC).  Results  are  from  the  simulations  on  the  48  x  50  x  48  quadratic  NURBS  mesh  described 
earlier. 
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(a)  streamwise  (downwind) 


(b)  spanwise  (crosswind) 
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Fig.  5  Crosswind-vertical  variation  of  velocity  fluctuations  v/  (defined  in  Eq.  (26))  in  flow  with 
LC.  Results  are  from  the  simulation  on  the  48  x  50  x  48  quadratic  NURBS  mesh. 


As  seen  in  Fig.  5,  a  change  in  sign  of  surface  intensified  V3  (in 
panel  ( b )  generates  the  surface  convergence  of  the  cell,  which  in 
turns  leads  to  the  downwelling  limb  of  the  cell.  The  downwelling 
limb  is  the  full-depth  region  characterized  by  negative  v'2  (in  panel 
(c)).  This  region  is  depicted  in  the  sketch  shown  in  Fig.  1 .  Further¬ 
more,  the  upwelling  limb  (region  with  positive  v'2)  of  the  cell  is 
larger  in  crosswind  extent  than  the  downwelling  region  (region 
with  negative  v4)  in  agreement  with  the  field  measurement  of  full- 
depth  LC  in  Ref.  [16].  At  mid  depth  the  upwelling  limb  is  approx¬ 
imately  1.6  larger  than  the  down- welling  limb,  which  is  close  to 
the  1 .4  factor  measured  in  the  field.  The  downwelling  limb  coin¬ 
cides  with  a  region  of  bottom-and  surface-intensified  positive  v\ 
(in  panel  (a)).  Note  that  this  region  of  full-depth  positive  v(  leads 
to  an  enhanced  downwind  current  as  depicted  in  Fig.  1.  The 
enhancement  of  the  downwind  current  within  the  Langmuir  cell 
downwelling  region  is  by  a  factor  of  approximately  10  u*  near  the 
surface  and  near  the  bottom  of  the  water  column.  Finally,  the  one¬ 
cell  structure  in  the  flow  with  C-L  vortex  forcing  (Fig.  5)  is  signif¬ 


icantly  different  in  structure  and  magnitude  of  fluctuations  from 
the  two-cell  structure  obtained  in  the  flow  without  C-L  vortex 
forcing  (Fig.  6). 

In  Fig.  7,  the  instantaneous  velocity  fluctuations  have  been 
made  dimensional  with  the  wind  stress  friction  velocity  reported 
by  Gargett  and  Wells  in  Ref.  [16]  during  their  field  observations 
of  full-depth  Langmuir  cells.  Magnitudes  of  these  fluctuations  in 
our  LES  are  in  close  agreement  with  those  measured  in  the  field 
(shown  in  Fig.  8)  as  well  as  with  those  computed  using  the  spec¬ 
tral  method  of  Tejada-Martlnez  et  al.  [9,21],  In  both,  computa¬ 
tions  and  field  experiments,  instantaneous  streamwise  and 
spanwise  velocity  fluctuations  are  in  the  ±8  cm/s  range  and  the 
vertical  velocity  fluctuation  is  in  the  ±4  cm/s  range.  Note  that  the 
field  measurements  in  Ref.  [16]  were  made  using  a  bottom- 
mounted,  upward-facing  acoustic  Doppler  current  profiler 
(ADCP)  in  a  15 -meter  deep  water  column  off  the  southern  New 
Jersey  coast  undergoing  strong  wind  and  wave  forcing.  Mean 
wind  stress  was  0.1  N/m2  and  mean  wave  height  was  1  m.  The 
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(b)  spanwise  (crosswind) 
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Fig.  6  Crosswind-vertical  variation  of  partially  averaged  velocity  fluctuations  vj  (defined  in 
Eq.  (26))  in  flow  without  LC.  Results  are  from  the  simulation  on  the  48  x  50  x  48  quadratic 
NURBS  mesh. 
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Fig.  7  Crosswind-vertical  variation  of  velocity  instantaneous  velocity  fluctuations  uj  (in  cm/s). 
Results  are  from  the  simulation  on  the  48  x  50  x  48  quadratic  NURBS  mesh.  Computational 
velocities  have  been  made  dimensional  with  wind  stress  friction  velocity  recorded  in  the  field 
during  episodes  of  full-depth  LC  [6],  Field  measurements  were  made  in  a  15-meters  deep  water 
column  under  a  wind  stress  of  0.1  N/m2. 


ADCP  was  not  able  to  make  a  reliable  measurement  of  the  upper¬ 
most  15%  of  the  water  column.  Furthermore,  the  computations  do 
not  take  into  account  the  effect  of  wave  breaking  at  the  surface. 
Thus  comparison  between  field  measurements  and  the  LES  should 
not  include  the  near-surface  region.  Comparison  of  panels  ( b )  in 
Figs.  7  and  8  shows  that  the  LES  is  able  to  resolve  the  near¬ 
bottom  intensification  of  the  full-depth  region  of  positive  down¬ 
wind  velocity  fluctuations  measured  in  the  field.  Furthermore,  in 
Fig.  8  note  that  the  region  of  downwelling  (panel  (c))  coincides 
with  a  region  of  positive  downwind  velocity  fluctuations  (panel 
a),  which  as  described  earlier  is  also  the  case  in  the  LES  (see  Figs. 
5  and  7). 

In  conclusion,  predictions  from  our  computation  with  C-L  vor¬ 
tex  forcing  compare  favorably  with  field  measurements  in  Ref. 


[16]  in  spite  of  the  low  Reynolds  number  of  the  computation 
(Re  =  395)  compared  to  the  Reynolds  number  of  the  observations 
(Re  ss  50,000). 

4.2  Mesh  Convergence.  Next  we  present  mesh  convergence 
studies  on  quadratic  NURBS  meshes  in  terms  of  mean  downwind 
velocity  and  turbulent  kinetic  energy  (TKE)  for  flow  with  full- 
depth  LC.  Mean  velocity,  TKE,  budgets  of  TKE  and  budgets  of 
TKE  components  for  this  flow  and  the  corresponding  flow  without 
LC  have  been  analyzed  in  detail  in  Ref.  [9].  Here  we  focus  strictly 
on  mesh  convergence.  Details  of  the  meshes  considered  are  given 
in  Table  1.  Note  that  for  the  coarsest  mesh  of  24x26x24 
elements,  the  first  wall-normal  mesh  knot  is  at  a  distance  yf 
=  uzAy/v  =  4.62.  For  the  the  finest  mesh  of  64  x  66  x  64  elements, 
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Fig.  8  Crosswind-vertical  variation  of  velocity  instantaneous  velocity  fluctuations  u/  (in  cm/s) 
during  episode  of  full-depth  Langmuir  cells  measured  during  field  experiments  of  Gargett  and 
Wells  [6]  using  a  bottom-mounted,  upward-facing  acoustic  Doppler  current  profiler  (ADCP). 
Field  measurements  were  made  in  a  15  ms-deep  water  column  under  a  wind  stress  of  0.1  N/m. 
This  figure  is  courtesy  of  Ann  Gargett. 
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Fig.  9  Convergence  of  mean  downwind  velocity  in  flow  with  LC.  Quadratic  NURBS  meshes 
were  used  for  all  cases. 


y j"  =  1.55.  Thus  for  all  meshes  considered,  the  first  point  off  the 
bottom  (top)  is  within  the  viscous  sublayer. 

Mean  downwind  velocity  is  expressed  as  {u),  recalling  that 
brackets  denote  averaging  over  time  and  downwind  and  cross¬ 


wind  directions.  TKE  is  defined  in  terms  of  velocity  fluctuations 
as  TKE  =  (u'u1  +  vV+  w'w')/2,  where  velocity  fluctuations  are 
again  obtained  via  the  classical  Reynolds  decomposition: 
u  =  (a)  +  u' . 


Fig.  10  Convergence  of  turbulent  kinetic  energy  (TKE)  in  flow  with  LC.  Quadratic  NURBS 
meshes  were  used  for  all  cases. 
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Fig.  11  Mean  downwind  velocity  in  flows  with  and  without  LC.  The  48x  50x48  quadratic 
NURBS  mesh  was  used  for  both  flows. 


Fig.  12  Mean  downwind  velocity  versus  wall-normal  (vertical)  direction  in  wall  (plus)  units  in 
flows  with  and  without  LC.  The  48  x  50  x  48  quadratic  NURBS  mesh  was  used  for  both  flows. 
Note  that  y+  =  uTylv. 
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We  only  present  our  computational  results,  because  no  direct 
numerical  simulation  (DNS)  of  this  test  case  exists  that  we  could 
use  as  a  benchmark  solution.  Figure  9  shows  convergence  of  the 
mean  velocity  profile.  The  24  x  26  x  24  mesh  gives  a  significant 
over-prediction  of  the  mean  flow.  The  results  improve  for  the  32 
x  34  x  32  mesh.  Further  improvement  is  seen  for  the  48  x  50 
x  48  mesh.  The  64  x  66  x  64  mesh  yields  a  nearly  indistinguish¬ 
able  mean  velocity  profile  from  the  48  x  50  x  48  case.  Similar 
convergence  pattern  is  observed  for  the  TKE  in  Fig.  10,  however, 
very  small  differences  between  the  48  x  50  x  48  and  the  64  x  66 
x  64  cases  are  visible  in  the  figure. 

4.3  Disruption  of  the  Log  Layer.  Figures  11  and  12  provide 
a  comparison  between  the  flow  with  LC  and  the  same  flow  with¬ 
out  LC  in  terms  of  mean  velocity  in  order  to  highlight  the  effects 
induced  by  LC.  The  action  of  LC  serves  to  homogenize  momen¬ 
tum  throughout  the  water  column  leading  to  a  near  constant  veloc¬ 
ity  profile  in  the  core  region  and  thinner  viscous  sublayers  at  the 
surface  and  bottom.  Figure  12  shows  mean  velocity  versus  wall- 
normal  direction  in  wall  units  in  the  lower  half  of  the  water  col¬ 
umn.  The  flow  without  LC  exhibits  a  well-developed,  near-bottom 
log-layer.  Meanwhile  in  the  flow  with  LC,  enhanced  mixing  asso¬ 
ciated  with  the  Langmuir  cells  disrupts  the  classical  log-layer 
region  inducing  an  extended  wake  region  at  depths  normally  char¬ 
acterized  by  the  log-law.  A  similiar  log-layer  disruption  has  been 
reported  in  Ref.  [21]  in  their  LES  of  full-depth  LC.  Disruption  of 
the  log-layer  by  the  action  of  LC  has  important  implications  for 
general  coastal  ocean  circulation  models  (GCOCMs).  Traditional 
Reynolds  Navier-Stokes  (RANS)  parameterizations  of  the  turbu¬ 
lent  bottom  boundary  layer  in  GCOCMs  assume  the  presence  of  a 
well-developed  log-layer.  Thus,  these  parameterizations  are  not 
able  to  properly  account  for  log-layer  disruption  caused  by  full- 
depth  LC  and  ultimately  wave-current  interaction. 

5  Conclusions  and  Future  Work 

We  have  successfully  extended  variational  multiscale  turbu¬ 
lence  modeling  procedures  to  the  C-L  equations  using  isogeomet¬ 
ric  analysis  based  on  quadratic  NURBS.  The  C-L  equations  were 
expressed  in  semilinear  form  revealing  an  advection-diffusion 
system  characterized  by  nonsymmetric  advective  matrices.  The 
weak  form  of  this  system  was  treated  with  the  residual-based  mul¬ 
tiscale  formulation  in  Ref.  [20]  together  with  stabilization  parame¬ 
ters  defined  in  terms  of  the  aforementioned  advective  matrices 
based  on  the  theory  presented  in  Ref.  [29].  No-slip  conditions 
were  enforced  weakly  following  the  approach  described  in  Refs. 
[23,35,36]. 

The  methodology  showed  good  convergence  properties  for  a 
wind-driven  shear  flow  characterized  by  large-scale  stream-wise 
structures  (full-depth  Langmuir  circulation)  in  agreement  with  the 
spectral  results  in  Ref.  [9]  and  the  field  measurements  in  Refs. 
[15,16].  A  major  impact  of  the  full-depth  Langmuir  cells  was 
shown  to  be  enhanced  mixing  of  momentum  leading  to  a  disrup¬ 
tion  of  the  classical  near-bottom  log-layer. 

Future  research  will  use  the  method  developed  in  this  manu¬ 
script  to  investigate  effects  of  complex  topography  on  Langmuir 
circulation  and  other  mixing  processes  in  coastal  oceans  and 
estuaries.  Furthermore,  we  will  explore  weak  enforcement  of 
Dirichlet  boundary  conditions  potentially  serving  as  a  model  of 
prohibitively  expensive  viscous  sublayers  in  continental  shelf 
flows. 
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